function udt = h2eom(t, u, params)
%   choice of Pi: 
%   Pi_xx =  2 P1 (1+UX^2) +  2 UX UY P2
%   Pi_yy = -2 P1 (1+UY^2) +  2 UX UY P2
%   Pi_xy =  P2 ( 2 + UX^2 + UY^2 )

    A = params.A;
    p1 = A+1;
    p2 = 2.*A;
    Nx = params.Nx;
    Ny = params.Ny;
    Nt = Nx.*Ny;

    %u = gpuArray(uu);
    


    % unpack variables 
    T  = reshape(u(0*Nt+1 : 1*Nt), Nx, Ny);
    ux = reshape(u(1*Nt+1 : 2*Nt), Nx, Ny);
    uy = reshape(u(2*Nt+1 : 3*Nt), Nx, Ny);
    P1 = reshape(u(3*Nt+1 : 4*Nt), Nx, Ny);
    P2 = reshape(u(4*Nt+1 : 5*Nt), Nx, Ny);

    ut = sqrt(1+ux.^2+uy.^2);
    

    global timer1 timer2 timer3 timer4 timer5;

    tic;

    Tdx  = dx(T);
    Tdy  = dy(T);
    uxdx = dx(ux);
    uxdy = dy(ux);
    uydx = dx(uy);
    uydy = dy(uy);
    P1dx = dx(P1);
    P1dy = dy(P1);
    P2dx = dx(P2);
    P2dy = dy(P2);

    utdx = (1/2).*ut.^(-1).*(2.*ux.*uxdx+2.*uy.*uydx);
    utdy = (1/2).*ut.^(-1).*(2.*ux.*uxdy+2.*uy.*uydy);


    ut2 = ut.^2;
    ut3 = ut.^3;
    ut4 = ut.^4;
    ux2 = ux.^2;
    ux3 = ux.^3;
    ux4 = ux.^4;
    uy2 = uy.^2;
    uy3 = uy.^3;
    uy4 = uy.^4;
    T2 = T.^2;
    T3 = T.^3;
    T4 = T.^4;

    timer1 = timer1  + toc();
    tic;

    % before ut simpl:
    % m11 = (-3).*T.^2.*(1+ux.^2+uy.^2).^(1/2);
    % m12 = (1/2).*(1+ux.^2+uy.^2).^(-1/2).*((-3).*T.^3.*ux+(-4).*P1.*ux.*(1+uy.^2)+(-2).*P2.*uy.*(2+(-1).*ux.^2+uy.^2));
    % m13 = (1/2).*(1+ux.^2+uy.^2).^(-1/2).*(((-3).*T.^3+4.*P1.*(1+ux.^2)).*uy+(-2).*P2.*ux.*(2+ux.^2+(-1).*uy.^2));

    % m22 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(2.*P2.*ux.*uy.*(2+3.*ux.^2+3.*ux.^4+uy.^2+(-3).*uy.^4)+3.*T.^3.*(ux.^4+3.*ux.^2.*(1+uy.^2)+2.*(1+uy.^2).^2)+4.*P1.*(ux.^4.*(1+(-3).*uy.^2)+2.*(1+uy.^2)+(-3).*ux.^2.*((-1)+uy.^4)));
    % m23 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(ux.*uy.*((-3).*T.^3.*(1+ux.^2+uy.^2)+4.*P1.*(1+ux.^2).*(1+3.*ux.^2+3.*uy.^2))+P2.*(8+(-6).*ux.^4+(-6).*ux.^6+12.*uy.^2+8.*uy.^4+ux.^2.*(8+6.*uy.^2+6.*uy.^4)));
    % m24 = 2.*(ux+ux.^3).*(1+ux.^2+uy.^2).^(-1/2);
    % m25 = uy.*(1+ux.^2+uy.^2).^(-1/2).*(2+3.*ux.^2+uy.^2);

    % m32 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(2.*P2.*(4+4.*uy.^2+(-3).*uy.^4+(-3).*uy.^6+3.*ux.^2.*(2+uy.^2)+ux.^4.*(4+3.*uy.^2))+(-1).*ux.*uy.*(3.*T.^3.*(1+ux.^2+uy.^2)+4.*P1.*(1+uy.^2).*(1+3.*ux.^2+3.*uy.^2)));
    % m33 = (1/4).*(1+ux.^2+uy.^2).^(-3/2).*(2.*P2.*ux.*uy.*(2+ux.^2+(-3).*ux.^4+3.*uy.^2+3.*uy.^4)+3.*T.^3.*(2+2.*ux.^4+3.*uy.^2+uy.^4+ux.^2.*(4+3.*uy.^2))+4.*P1.*((-2)+(-3).*uy.^2+3.*ux.^4.*uy.^2+(-1).*uy.^4+ux.^2.*((-2)+3.*uy.^4)));
    % m34 = (-2).*(1+ux.^2+uy.^2).^(-1/2).*(uy+uy.^3);
    % m35 = ux.*(1+ux.^2+uy.^2).^(-1/2).*(2+ux.^2+3.*uy.^2);

    % m42 = (1/12).*T.^(-4).*(1+ux.^2+uy.^2).^(-1/2).*(16.*(1+A).*P1.^2.*ux.*(1+uy.^2).*(2+ux.^2+uy.^2)+3.*T.^3.*(T.^3.*ux.*(2+ux.^2+3.*uy.^2)+2.*P2.*uy.*((2+A).*ux.^4+(-2).*ux.^2.*((-1)+A+2.*uy.^2+A.*uy.^2)+(2+uy.^2).*(2+(2+A).*uy.^2)))+2.*P1.*(4.*(1+A).*P2.*uy.*((-1).*ux.^4+(2+uy.^2).^2)+3.*T.^3.*ux.*(6+7.*uy.^2+4.*uy.^4+ux.^2.*(3+(-4).*uy.^2)+A.*(6+7.*uy.^2+2.*uy.^4+ux.^2.*(3+(-2).*uy.^2)))));
    % m43 = (1/12).*T.^(-4).*(1+ux.^2+uy.^2).^(-1/2).*((-16).*(1+A).*P1.^2.*(1+ux.^2).*uy.*(2+ux.^2+uy.^2)+(-3).*T.^3.*(T.^3.*uy.*(2+3.*ux.^2+uy.^2)+2.*P2.*ux.*(4+(2+A).*ux.^4+(-2).*((-1)+A).*uy.^2+(2+A).*uy.^4+ux.^2.*(6+(-4).*uy.^2+(-2).*A.*((-1)+uy.^2))))+2.*P1.*(4.*(1+A).*P2.*ux.*(4+4.*ux.^2+ux.^4+(-1).*uy.^4)+3.*T.^3.*uy.*(4.*ux.^4+ux.^2.*(7+(-4).*uy.^2)+3.*(2+uy.^2)+A.*(2.*ux.^4+ux.^2.*(7+(-2).*uy.^2)+3.*(2+uy.^2)))));
    % m44 = (1+A).*T.^(-1).*(1+ux.^2+uy.^2).^(1/2).*(2+ux.^2+uy.^2);

    % m52 = (1/6).*T.^(-4).*(1+ux.^2+uy.^2).^(-1/2).*(3.*T.^6.*uy.*(1+uy.^2)+4.*(1+A).*P2.^2.*uy.*((-1).*ux.^4+(2+uy.^2).^2)+4.*P1.*(1+uy.^2).*((-3).*T.^3.*(1+(2+A).*ux.^2).*uy+2.*(1+A).*P2.*ux.*(2+ux.^2+uy.^2))+3.*P2.*T.^3.*ux.*(6+uy.^2+(-4).*uy.^4+ux.^2.*(5+4.*uy.^2)+A.*(6+uy.^2+(-2).*uy.^4+ux.^2.*(5+2.*uy.^2))));
    % m53 = (1/6).*T.^(-4).*(1+ux.^2+uy.^2).^(-1/2).*(4.*(1+A).*P2.^2.*ux.*(4+4.*ux.^2+ux.^4+(-1).*uy.^4)+3.*T.^3.*ux.*(1+ux.^2).*(T.^3+4.*P1.*(1+(2+A).*uy.^2))+(-1).*P2.*uy.*(8.*(1+A).*P1.*(1+ux.^2).*(2+ux.^2+uy.^2)+3.*T.^3.*((-6)+4.*ux.^4+(-5).*uy.^2+(-1).*ux.^2.*(1+4.*uy.^2)+A.*((-6)+2.*ux.^4+(-5).*uy.^2+(-1).*ux.^2.*(1+2.*uy.^2)))));
    % m55 = (1+A).*T.^(-1).*(1+ux.^2+uy.^2).^(1/2).*(2+ux.^2+uy.^2);
    
    % after ut simp:
    m11 = (-3).*T2.*ut;
    m12 = (-1/2).*ut.^(-1).*(3.*T3.*ux+4.*P1.*ux.*(1+uy2)+2.*P2.*uy.*(2+(-1).*ux2+uy2));
    m13 = (1/2).*ut.^(-1).*(((-3).*T3+4.*P1.*(1+ux2)).*uy+(-2).*P2.*ux.*(1+ut2+(-2).*uy2));

    m22 = (1/4).*ut.^(-3).*(3.*T3.*ut2.*(1+ut2+uy2)+4.*P1.*(2.*(1+uy2)+ux2.*(2+ut2+(2+(-3).*ut2).*uy2))+2.*P2.*ux.*uy.*(2+3.*(ux2+ux4)+uy2+(-3).*uy4));
    m23 = (1/4).*ut.^(-3).*(ux.*uy.*(ut2.*((-3).*T3+P1.*(4+12.*ux2))+8.*P1.*uy2)+P2.*(ut2.*(8+(-6).*ux4+6.*ux2.*uy2)+4.*(uy2+2.*uy4)));
    m24 = 2.*ut.^(-1).*(ux+ux3);
    m25 = ut.^(-1).*uy.*(2+3.*ux2+uy2);

    m32 = (1/4).*ut.^(-3).*(ux.*uy.*((-3).*T3.*ut2+(-4).*P1.*((-2)+3.*ut2).*(1+uy2))+P2.*(4+12.*uy2+ut4.*(8+6.*uy2)+ut2.*((-4)+(-22).*uy2+(-12).*uy4)+8.*uy4));
    m33 = (1/4).*ut.^(-3).*(3.*T3.*ut2.*(2+2.*ux2+uy2)+2.*P2.*ux.*uy.*(2+ux2+(-3).*ux4+3.*(uy2+uy4))+P1.*(ut2.*((-8)+12.*((-1)+ux2).*uy2)+8.*(uy2+uy4)));
    m34 = (-2).*ut.^(-1).*(uy+uy3);
    m35 = ut.^(-1).*ux.*(1+ut2+2.*uy2);

    m42 = (1/12).*T.^(-4).*ut.^(-1).*(3.*T3.*(T3.*ux.*(1+ut2+2.*uy2)+(-1).*p2.*uy.*((-2).*P1.*ux.*uy.*(ux2+(-1).*uy2)+P2.*(ux4+(-2).*ux2.*((-2)+uy2)+(2+uy2).^2)))+p1.*(16.*P1.^2.*(1+ut2).*ux.*(1+uy2)+12.*P2.*T3.*uy.*(2+ux2+ux4+3.*uy2+(-2).*ux2.*uy2+uy4)+P1.*(8.*P2.*uy.*((-1).*ux4+(2+uy2).^2)+(-6).*T3.*ux.*((-3)+ut2.*((-3)+4.*uy2)+(-8).*(uy2+uy4)))));
    m43 = (1/12).*T.^(-4).*ut.^(-1).*(3.*T3.*((-1).*T3.*uy.*(2+3.*ux2+uy2)+p2.*ux.*((-2).*P1.*ux.*uy.*(ux2+(-1).*uy2)+P2.*(ux4+(-2).*ux2.*((-2)+uy2)+(2+uy2).^2)))+2.*p1.*((-8).*P1.^2.*(1+ut2).*(1+ux2).*uy+P1.*(3.*T3.*uy.*(4.*ux4+ux2.*(7+(-4).*uy2)+3.*(2+uy2))+4.*P2.*ux.*((2+ux2).^2+(-1).*uy4))+(-6).*P2.*T3.*ux.*(2+ux4+ux2.*(3+(-2).*uy2)+uy2+uy4)));
    m44 = p1.*T.^(-1).*(ut+ut3);

    m52 = (1/6).*T.^(-4).*ut.^(-1).*(3.*T3.*uy.*(T3.*(1+uy2)+2.*P1.*p2.*(1+ux2).*(1+uy2)+p2.*P2.*ux.*uy.*((-1).*ux2+uy2))+p1.*(4.*P1.*(2.*P2.*(1+ut2).*ux+(-3).*T3.*(1+2.*ux2).*uy).*(1+uy2)+P2.*(4.*P2.*uy.*((-1).*ux4+(2+uy2).^2)+3.*T3.*ux.*(6+uy2+ux2.*(5+4.*uy2)+(-4).*uy4))));
    m53 = (1/6).*T.^(-4).*ut.^(-1).*(3.*T3.*ux.*(T3.*(1+ux2)+p2.*P2.*ux.*uy.*(ux2+(-1).*uy2)+(-2).*P1.*p2.*(1+ux2).*(1+uy2))+p1.*(12.*P1.*T3.*(ux+ux3).*(1+2.*uy2)+P2.*uy.*((-8).*P1.*(1+ut2).*(1+ux2)+3.*T3.*(6+ux2+(-4).*ux4+(5+4.*ux2).*uy2))+4.*P2.^2.*ux.*((2+ux2).^2+(-1).*uy4)));
    m55 = p1.*T.^(-1).*(ut+ut3);




    timer2 = timer2 + toc;
    tic;



    % before round ut simplification:
    % R1 = (-3).*T.^2.*Tdx.*ux+(-3).*T.^2.*Tdy.*uy+(1/2).*uxdx.*(1+ux.^2+uy.^2).^(-1).*(2.*P2.*ux.*(ux+(-1).*uy).*uy.*(ux+uy)+(-4).*P1.*(1+ux.^2).*(1+uy.^2)+(-3).*T.^3.*(1+ux.^2+uy.^2))+(-1).*uxdy.*(1+ux.^2+uy.^2).^(-1).*(2.*P1.*ux.*uy.*(1+uy.^2)+P2.*(2+ux.^2+(-1).*((-3)+ux.^2).*uy.^2+uy.^4))+(1+ux.^2+uy.^2).^(-1).*(2.*P1.*ux.*(1+ux.^2).*uy+(-1).*P2.*(2+ux.^4+uy.^2+(-1).*ux.^2.*((-3)+uy.^2))).*uydx+(1/2).*(1+ux.^2+uy.^2).^(-1).*(4.*P1.*(1+ux.^2).*(1+uy.^2)+2.*P2.*ux.*uy.*((-1).*ux.^2+uy.^2)+(-3).*T.^3.*(1+ux.^2+uy.^2)).*uydy;
    % R2 = (3/2).*T.^2.*Tdx+2.*P1dx.*(1+ux.^2)+2.*P2dx.*ux.*uy+P2dy.*(2+ux.^2+uy.^2)+(1/4).*uxdx.*(1+ux.^2+uy.^2).^(-1).*(3.*T.^3.*ux.*(1+ux.^2+uy.^2)+4.*P1.*ux.*(1+uy.^2+ux.^2.*(1+(-3).*uy.^2))+2.*P2.*uy.*(ux.^2.*(4+3.*ux.^2+(-3).*uy.^2)+4.*(1+uy.^2)))+(1/2).*uxdy.*(1+ux.^2+uy.^2).^(-1).*(P2.*ux.*((-2)+ux.^2+((-5)+3.*ux.^2).*uy.^2+(-3).*uy.^4)+3.*uy.*((-2).*P1.*ux.^2.*(1+uy.^2)+T.^3.*(1+ux.^2+uy.^2)))+(1/2).*ux.*(1+ux.^2+uy.^2).^(-1).*(6.*P1.*ux.*(1+ux.^2).*uy+P2.*((-2)+(-3).*ux.^4+uy.^2+ux.^2.*((-5)+3.*uy.^2))).*uydx+(1/4).*(1+ux.^2+uy.^2).^(-1).*(12.*P1.*(ux+ux.^3).*(1+uy.^2)+(-3).*T.^3.*ux.*(1+ux.^2+uy.^2)+2.*P2.*uy.*(4+4.*ux.^2+(-3).*ux.^4+(4+3.*ux.^2).*uy.^2)).*uydy;
    % R3 = (3/2).*T.^2.*Tdy+2.*P2dy.*ux.*uy+(-2).*P1dy.*(1+uy.^2)+P2dx.*(2+ux.^2+uy.^2)+(1/2).*uxdy.*uy.*(1+ux.^2+uy.^2).^(-1).*((-6).*P1.*ux.*uy.*(1+uy.^2)+P2.*((-2)+ux.^2+((-5)+3.*ux.^2).*uy.^2+(-3).*uy.^4))+(1/4).*uxdx.*(1+ux.^2+uy.^2).^(-1).*((-3).*uy.*(4.*P1.*(1+ux.^2).*(1+uy.^2)+T.^3.*(1+ux.^2+uy.^2))+2.*P2.*ux.*(4+4.*uy.^2+(-3).*uy.^4+ux.^2.*(4+3.*uy.^2)))+(1/2).*(1+ux.^2+uy.^2).^(-1).*(3.*T.^3.*ux.*(1+ux.^2+uy.^2)+uy.*(6.*P1.*ux.*(1+ux.^2).*uy+P2.*((-2)+uy.^2+ux.^2.*((-5)+(-3).*ux.^2+3.*uy.^2)))).*uydx+(1/4).*(1+ux.^2+uy.^2).^(-1).*((-4).*P1.*(1+ux.^2).*uy+4.*P1.*((-1)+3.*ux.^2).*uy.^3+3.*T.^3.*uy.*(1+ux.^2+uy.^2)+2.*P2.*ux.*(4+4.*uy.^2+3.*uy.^4+ux.^2.*(4+(-3).*uy.^2))).*uydy;
    % R4 = P1.*(2+ux.^2+uy.^2)+(1+A).*P1dx.*T.^(-1).*ux.*(2+ux.^2+uy.^2)+(1+A).*P1dy.*T.^(-1).*uy.*(2+ux.^2+uy.^2)+(1/6).*T.^(-4).*uxdy.*(1+ux.^2+uy.^2).^(-1).*(8.*(1+A).*P1.^2.*ux.*uy.*(1+uy.^2).*(2+ux.^2+uy.^2)+3.*T.^3.*((-1).*A.*P2.*(2+ux.^2).^2+T.^3.*ux.*(1+ux.^2).*uy+P2.*(4+(-4).*A+2.*ux.^2+(2+A).*ux.^4).*uy.^2+T.^3.*ux.*uy.^3+P2.*(6+A+(-2).*(2+A).*ux.^2).*uy.^4+(2+A).*P2.*uy.^6)+2.*P1.*(3.*T.^3.*uy.*(A.*ux.^3+(-1).*(2+A).*ux.*((-1)+ux.^2).*uy.^2+(2+A).*ux.*uy.^4)+2.*(1+A).*P2.*((2+ux.^2).^2+(8+2.*ux.^2+(-1).*ux.^4).*uy.^2+5.*uy.^4+uy.^6)))+(1/12).*T.^(-4).*uxdx.*(1+ux.^2+uy.^2).^(-1).*(16.*(1+A).*P1.^2.*(1+ux.^2).*(1+uy.^2).*(2+ux.^2+uy.^2)+2.*P1.*((-4).*(1+A).*P2.*ux.*(ux+(-1).*uy).*uy.*(ux+uy).*(2+ux.^2+uy.^2)+3.*T.^3.*(3.*(1+A).*(2+3.*ux.^2+ux.^4)+(9.*(1+A)+10.*(1+A).*ux.^2+(-2).*(2+A).*ux.^4).*uy.^2+(3.*(1+A)+2.*(2+A).*ux.^2).*uy.^4))+3.*T.^3.*(T.^3.*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2)+2.*P2.*ux.*uy.*((2+A).*ux.^4+(2+uy.^2).*(2+(2+A).*uy.^2)+(-2).*ux.^2.*((-1)+A+(2+A).*uy.^2))))+(1/6).*T.^(-4).*(1+ux.^2+uy.^2).^(-1).*((-8).*(1+A).*P1.^2.*ux.*(1+ux.^2).*uy.*(2+ux.^2+uy.^2)+2.*P1.*(2.*(1+A).*P2.*(2+ux.^2+uy.^2).*(2+ux.^4+uy.^2+(-1).*ux.^2.*((-3)+uy.^2))+3.*T.^3.*ux.*uy.*((2+A).*ux.^4+A.*uy.^2+(-1).*(2+A).*ux.^2.*((-1)+uy.^2)))+(-3).*T.^3.*(A.*P2.*((-2)+ux.^2+(-1).*uy.^2).*(2+ux.^4+uy.^2+(-1).*ux.^2.*((-3)+uy.^2))+ux.*(T.^3.*uy.*(1+ux.^2+uy.^2)+2.*P2.*ux.*(2+ux.^4+uy.^2+uy.^4+ux.^2.*(3+(-2).*uy.^2))))).*uydx+(1/12).*T.^(-4).*(1+ux.^2+uy.^2).^(-1).*((-16).*(1+A).*P1.^2.*(1+ux.^2).*(1+uy.^2).*(2+ux.^2+uy.^2)+(-3).*T.^3.*(T.^3.*(1+ux.^2+uy.^2).*(2+ux.^2+uy.^2)+2.*P2.*ux.*uy.*(4+2.*(3+A).*ux.^2+(2+A).*ux.^4+(-2).*((-1)+A+(2+A).*ux.^2).*uy.^2+(2+A).*uy.^4))+2.*P1.*(4.*(1+A).*P2.*ux.*(ux+(-1).*uy).*uy.*(ux+uy).*(2+ux.^2+uy.^2)+3.*T.^3.*(3.*(1+A).*(2+3.*ux.^2+ux.^4)+(9.*(1+A)+10.*(1+A).*ux.^2+2.*(2+A).*ux.^4).*uy.^2+(3.*(1+A)+(-2).*(2+A).*ux.^2).*uy.^4))).*uydy;
    % R5 = P2.*(2+ux.^2+uy.^2)+(1+A).*P2dx.*T.^(-1).*ux.*(2+ux.^2+uy.^2)+(1+A).*P2dy.*T.^(-1).*uy.*(2+ux.^2+uy.^2)+(1/6).*T.^(-4).*uxdy.*(1+ux.^2+uy.^2).^(-1).*(3.*T.^3.*(1+uy.^2).*((4.*A.*P1+T.^3).*(1+ux.^2)+((-4).*P1+T.^3+(-4).*(2+A).*P1.*ux.^2).*uy.^2)+(-4).*(1+A).*P2.^2.*(2+ux.^2+uy.^2).*(ux.^2.*((-1)+uy.^2)+(-1).*(1+uy.^2).*(2+uy.^2))+2.*P2.*ux.*uy.*(4.*(1+A).*P1.*(1+uy.^2).*(2+ux.^2+uy.^2)+3.*T.^3.*(ux+(-1).*uy).*(ux+uy).*(1+(2+A).*uy.^2)))+(1/6).*T.^(-4).*uxdx.*(1+ux.^2+uy.^2).^(-1).*(4.*P1.*(1+uy.^2).*((-3).*T.^3.*ux.*(1+(2+A).*ux.^2).*uy+2.*(1+A).*P2.*(1+ux.^2).*(2+ux.^2+uy.^2))+P2.*((-4).*(1+A).*P2.*ux.*(ux+(-1).*uy).*uy.*(ux+uy).*(2+ux.^2+uy.^2)+3.*T.^3.*((1+A).*(6+9.*ux.^2+5.*ux.^4)+(9.*(1+A)+4.*(1+A).*ux.^2+2.*(2+A).*ux.^4).*uy.^2+(3.*(1+A)+(-2).*(2+A).*ux.^2).*uy.^4)))+(1/6).*T.^(-4).*(1+ux.^2+uy.^2).^(-1).*(4.*(1+A).*P2.^2.*(2+ux.^2+uy.^2).*(2+ux.^4+uy.^2+(-1).*ux.^2.*((-3)+uy.^2))+(-2).*P2.*ux.*uy.*(3.*T.^3.*(1+(2+A).*ux.^2).*(ux+(-1).*uy).*(ux+uy)+4.*(1+A).*P1.*(1+ux.^2).*(2+ux.^2+uy.^2))+3.*T.^3.*(1+ux.^2).*(T.^3.*(1+ux.^2+uy.^2)+4.*P1.*ux.^2.*(1+2.*uy.^2)+4.*A.*P1.*((-1)+((-1)+ux.^2).*uy.^2))).*uydx+(1/6).*T.^(-4).*(1+ux.^2+uy.^2).^(-1).*(4.*P1.*(1+ux.^2).*((-2).*(1+A).*P2.*(1+uy.^2).*(2+ux.^2+uy.^2)+3.*T.^3.*ux.*uy.*(1+(2+A).*uy.^2))+P2.*(4.*(1+A).*P2.*ux.*(ux+(-1).*uy).*uy.*(ux+uy).*(2+ux.^2+uy.^2)+3.*T.^3.*(3.*(1+A).*(2+3.*ux.^2+ux.^4)+(9.*(1+A)+4.*(1+A).*ux.^2+(-2).*(2+A).*ux.^4).*uy.^2+(5.*(1+A)+2.*(2+A).*ux.^2).*uy.^4))).*uydy;


    R1 = (-3).*T2.*(Tdx.*ux+Tdy.*uy)+(-3/2).*T3.*((-1).*((-1)+ux2).*uxdx+ut.*(utdx.*ux+utdy.*uy)+uydy+(-1).*uy.*(ux.*(uxdy+uydx)+uy.*uydy))+(-2).*P1.*ut.^(-1).*((-1).*utdx.*(ux+ux3)+utdy.*(uy+uy3)+ut.*((1+ux2).*uxdx+(-1).*(1+uy2).*uydy))+P2.*ut.^(-1).*(utdx.*uy.*(2+3.*ux2+uy2)+utdy.*ux.*(1+ut2+2.*uy2)+(-1).*ut3.*(uxdy+uydx)+(-1).*ut.*(uxdy+uydx+2.*ux.*uy.*(uxdx+uydy)));
    R2 = (3/2).*T2.*Tdx+P2dy.*(1+ut2)+2.*(P1dx+P1dx.*ux2+P2dx.*ux.*uy)+(1/2).*P2.*ut.^(-1).*(3.*utdx.*ux.*uy.*(2+3.*ux2+uy2)+utdy.*(3.*ux4+8.*(1+uy2)+ux2.*(14+9.*uy2))+(-3).*ut3.*ux.*(uxdy+uydx)+ut.*(ux.*((-7).*uxdy+uydx)+4.*uy.*(uxdx+(-1).*uydy)+(-6).*ux2.*uy.*(uxdx+uydy)))+(-3/4).*T3.*((-2).*uxdy.*uy+3.*ut.*ux.*(utdx.*ux+utdy.*uy)+ux.*(((-1)+(-3).*ux2).*uxdx+uydy+(-3).*uy.*(ux.*(uxdy+uydx)+uy.*uydy)))+P1.*ut.^(-1).*ux.^(-1).*uy.^(-1).*(3.*uy.*(utdx.*(ux.^5+ux3)+utdy.*uy.*(1+uy2).^2)+ut2.*utdy.*(2.*ux2+(-5).*uy2+(-3).*uy4)+ut.*((-2).*ux3.*uxdy+(-3).*ux4.*uxdx.*uy+2.*ux.*uxdy.*uy2+2.*uy3.*uydy+ux2.*uy.*(uxdx+uydy+3.*uy2.*uydy)));
    R3 = P2dx+(3/2).*T2.*Tdy+P2dx.*ut2+2.*P2dy.*ux.*uy+(-2).*P1dy.*(1+uy2)+(-3/4).*T3.*(3.*ut.*uy.*(utdx.*ux+utdy.*uy)+uxdx.*(uy+(-3).*ux2.*uy)+ux.*((-2).*uydx+(-3).*uy2.*(uxdy+uydx))+uy.*((-1)+(-3).*uy2).*uydy)+P1.*ut.^(-1).*ux.^(-1).*uy.^(-1).*((-3).*utdy.*ux.*(1+uy2).*uy3+utdx.*(2.*(ux2+ux4)+((-2)+3.*(ux2+ux4)).*uy2+(-2).*uy4)+(-1).*ut.*(ux3.*uxdx.*(2+3.*uy2)+2.*ux2.*uy.*uydx+(-2).*uy3.*uydx+ux.*uy2.*(uxdx+uydy+(-3).*uy2.*uydy)))+(1/2).*P2.*ut.^(-1).*(ut2.*(3.*utdy.*ux.*uy+utdx.*(8+9.*uy2))+3.*(utdy.*ux+(-1).*utdx.*uy).*(uy+2.*uy3)+(-3).*ut3.*uy.*(uxdy+uydx)+ut.*(uy.*(uxdy+(-7).*uydx)+ux.*(uxdx.*((-4)+(-6).*uy2)+(-2).*((-2)+3.*uy2).*uydy)));
    R4 = P1.*(1+ut2)+(1/4).*p2.*T.^(-1).*ut.^(-1).*((-2).*P1.*ux.*uy.*(ux2+(-1).*uy2)+P2.*(ux4+(-2).*ux2.*((-2)+uy2)+(2+uy2).^2)).*(utdy.*ux+(-1).*utdx.*uy+ut.*((-1).*uxdy+uydx))+(1/4).*T2.*((-1).*ut3.*utdx.*ux+uxdx.*(2+ux2.*(3+ux2+(-1).*uy2)+uy2)+ut.*(utdy.*uy.*(2+(-1).*ux2+uy2)+utdx.*ux.*((-1)+2.*uy2))+ux.*uy.*(ux2+(-1).*uy2).*(uxdy+uydx)+ut2.*((-1)+uy2).*uydy+((-1)+(-3).*uy2+(-2).*uy4).*uydy)+p1.*((2/3).*P1.*T.^(-4).*ut.^(-1).*(P2.*(1+ut2).*((-1).*utdx.*uy.*(2+3.*ux2+uy2)+(-1).*utdy.*ux.*(1+ut2+2.*uy2)+ut3.*(uxdy+uydx)+ut.*(uxdy+uydx+2.*ux.*uy.*(uxdx+uydy)))+2.*P1.*(uy.*(utdy+utdx.*ux.*uy+utdy.*uy2)+ut2.*((-1).*utdx.*ux.*(2+ux2)+utdy.*(uy+uy3))+ut3.*((2+ux2).*uxdx+(-1).*(1+uy2).*uydy)+(-1).*ut.*(uydy+uy2.*(uxdx+uydy))))+T.^(-1).*((1+ut2).*(P1dx.*ux+P1dy.*uy)+P2.*ut.^(-1).*(2.*utdx.*ux2.*uy.*(ux2+(-1).*uy2)+2.*utdy.*ux.*(ux2+(-1).*uy2).*uy2+(-1).*ut3.*ux.*uy.*uydy+ut.*((-1).*ux3.*uxdx.*uy+uxdy.*uy2.*(2+uy2)+(-1).*ux4.*uydx+ux2.*((-1).*uxdy.*uy2+((-2)+uy2).*uydx)+ux.*uy.*(uxdx.*(2+uy2)+((-1)+2.*uy2).*uydy)))+(1/2).*P1.*ut.^(-1).*(2.*ut2.*(utdx.*ux+utdy.*uy).*(2+3.*ux2+(-1).*uy2)+8.*(utdx.*ux+utdy.*uy).*uy2.*(1+uy2)+ut3.*((-6).*ux.*uxdy.*uy+(3+(-2).*uy2).*uydy)+(-1).*ut.*(uxdx.*((-6)+ux2+6.*ux4+((-3)+2.*ux2).*uy2)+(-3).*uydy+2.*uy.*(ux.*(uxdy.*((-1)+(-2).*uy2)+(2+ux2+3.*uy2).*uydx)+(uy+2.*uy3).*uydy)))));
    R5 = P2.*(1+ut2)+(1/2).*p2.*T.^(-1).*ut.^(-1).*(2.*P1.*(1+ux2).*(1+uy2)+P2.*ux.*uy.*((-1).*ux2+uy2)).*((-1).*utdy.*ux+utdx.*uy+ut.*(uxdy+(-1).*uydx))+(1/2).*T2.*((-1).*ut.*(utdx.*(1+ux2).*uy+utdy.*ux.*(1+uy2))+(1+ux2).*(uxdy+ux.*uxdx.*uy+uxdy.*uy2+uydx+uy2.*uydx)+ux.*(uy+uy3).*uydy)+p1.*(T.^(-1).*((1+ut2).*(P2dx.*ux+P2dy.*uy)+(1/2).*P2.*ut.^(-1).*((-1).*ut.*((-11)+ut2+4.*ut4).*uxdx+12.*utdy.*uy+2.*(utdx.*ux.*(6+9.*ux2+2.*ux4+(9+8.*ux2).*uy2+2.*uy4)+uy.*(utdy.*(ux2.*(9+2.*ux2)+(9+8.*ux2).*uy2+2.*uy4)+2.*ut.*(uxdx.*uy.*(3+uy2)+ux.*((-1).*uxdy.*(2+ut2+uy2)+((-1)+(-2).*ut2+uy2).*uydx))))+ut.*(3+ut2.*(3+(-8).*uy2)+(-4).*uy2+4.*uy4).*uydy)+2.*P1.*ut.^(-1).*ux.^(-1).*uy.^(-1).*(utdx.*ux.*(1+ux2).*(ux2+(-1).*uy2).*(1+uy2)+utdy.*(1+ux2).*uy.*(ux2+(-1).*uy2).*(1+uy2)+ut.*((-1).*ux3.*(ux.*uxdx+uxdy.*uy).*(1+uy2)+ux.*(1+ux2).*uy3.*uydx+(1+ux2).*uy4.*uydy)))+T.^(-4).*((4/3).*P1.*P2.*(uxdx.*(ut2.*(2+ux2)+(-1).*uy2)+ut.^(-1).*(1+ut2).*utdy.*uy.*(1+uy2)+ut.^(-1).*utdx.*ux.*((-1).*ut2.*(2+ux2)+uy2)+(-1).*(1+ut2).*(1+uy2).*uydy)+(2/3).*P2.^2.*ut.^(-1).*(1+ut2).*((-1).*utdx.*uy.*(2+3.*ux2+uy2)+(-1).*utdy.*ux.*(1+ut2+2.*uy2)+ut3.*(uxdy+uydx)+ut.*(uxdy+uydx+2.*ux.*uy.*(uxdx+uydy)))));








    timer3 = timer3 + toc;
    tic;



    mdet = m11.*m25.*m34.*m43.*m52+(-1).*m11.*m24.*m35.*m43.*m52+(-1).*m11.*m25.*m33.*m44.*m52+m11.*m23.*m35.*m44.*m52+(-1).*m11.*m25.*m34.*m42.*m53+m11.*m24.*m35.*m42.*m53+m11.*m25.*m32.*m44.*m53+(-1).*m11.*m22.*m35.*m44.*m53+(-1).*m11.*m24.*m33.*m42.*m55+m11.*m23.*m34.*m42.*m55+m11.*m24.*m32.*m43.*m55+(-1).*m11.*m22.*m34.*m43.*m55+(-1).*m11.*m23.*m32.*m44.*m55+m11.*m22.*m33.*m44.*m55;

   
    i11 = m35.*((-1).*m24.*m43.*m52+m23.*m44.*m52+m24.*m42.*m53+(-1).*m22.*m44.*m53)+m25.*(m34.*m43.*m52+(-1).*m33.*m44.*m52+(-1).*m34.*m42.*m53+m32.*m44.*m53)+((-1).*m24.*m33.*m42+m23.*m34.*m42+m24.*m32.*m43+(-1).*m22.*m34.*m43+(-1).*m23.*m32.*m44+m22.*m33.*m44).*m55;
    i12 = (-1).*m13.*(m35.*m44.*m52+m34.*m42.*m55+(-1).*m32.*m44.*m55)+m12.*(m35.*m44.*m53+m34.*m43.*m55+(-1).*m33.*m44.*m55);
    i13 = m25.*m44.*(m13.*m52+(-1).*m12.*m53)+(m13.*m24.*m42+(-1).*m12.*m24.*m43+(-1).*m13.*m22.*m44+m12.*m23.*m44).*m55;
    i14 = (-1).*(m25.*m34+(-1).*m24.*m35).*(m13.*m52+(-1).*m12.*m53)+((-1).*m13.*m24.*m32+m12.*m24.*m33+m13.*m22.*m34+(-1).*m12.*m23.*m34).*m55;
    i15 = (m25.*m34+(-1).*m24.*m35).*(m13.*m42+(-1).*m12.*m43)+((-1).*m13.*m25.*m32+m12.*m25.*m33+m13.*m22.*m35+(-1).*m12.*m23.*m35).*m44;
    
    i22 = (-1).*m11.*(m35.*m44.*m53+m34.*m43.*m55+(-1).*m33.*m44.*m55);
    i23 = m11.*(m25.*m44.*m53+m24.*m43.*m55+(-1).*m23.*m44.*m55);
    i24 = m11.*((-1).*m25.*m34.*m53+m24.*m35.*m53+(-1).*m24.*m33.*m55+m23.*m34.*m55);
    i25 = m11.*(m25.*m34.*m43+(-1).*m24.*m35.*m43+(-1).*m25.*m33.*m44+m23.*m35.*m44);
    
    i32 = m11.*(m35.*m44.*m52+m34.*m42.*m55+(-1).*m32.*m44.*m55);
    i33 = (-1).*m11.*(m25.*m44.*m52+m24.*m42.*m55+(-1).*m22.*m44.*m55);
    i34 = m11.*(m25.*m34.*m52+(-1).*m24.*m35.*m52+m24.*m32.*m55+(-1).*m22.*m34.*m55);
    i35 = m11.*((-1).*m25.*m34.*m42+m24.*m35.*m42+m25.*m32.*m44+(-1).*m22.*m35.*m44);

    i42 = m11.*((-1).*m35.*m43.*m52+m35.*m42.*m53+(-1).*m33.*m42.*m55+m32.*m43.*m55);
    i43 = m11.*(m25.*m43.*m52+(-1).*m25.*m42.*m53+m23.*m42.*m55+(-1).*m22.*m43.*m55);
    i44 = m11.*((-1).*m25.*m33.*m52+m23.*m35.*m52+m25.*m32.*m53+(-1).*m22.*m35.*m53+(-1).*m23.*m32.*m55+m22.*m33.*m55);
    i45 = m11.*(m25.*m33.*m42+(-1).*m23.*m35.*m42+(-1).*m25.*m32.*m43+m22.*m35.*m43);

    i52 = m11.*(m34.*m43.*m52+(-1).*m33.*m44.*m52+(-1).*m34.*m42.*m53+m32.*m44.*m53);
    i53 = m11.*((-1).*m24.*m43.*m52+m23.*m44.*m52+m24.*m42.*m53+(-1).*m22.*m44.*m53);
    i54 = m11.*(m24.*m33.*m52+(-1).*m23.*m34.*m52+(-1).*m24.*m32.*m53+m22.*m34.*m53);
    i55 = m11.*((-1).*m24.*m33.*m42+m23.*m34.*m42+m24.*m32.*m43+(-1).*m22.*m34.*m43+(-1).*m23.*m32.*m44+m22.*m33.*m44);

    i11 = i11./mdet;
    i12 = i12./mdet;
    i13 = i13./mdet;
    i14 = i14./mdet;
    i15 = i15./mdet;

    i22 = i22./mdet;
    i23 = i23./mdet;
    i24 = i24./mdet;
    i25 = i25./mdet;

    i32 = i32./mdet;
    i33 = i33./mdet;
    i34 = i34./mdet;
    i35 = i35./mdet;

    i42 = i42./mdet;
    i43 = i43./mdet;
    i44 = i44./mdet;
    i45 = i45./mdet;

    i52 = i52./mdet;
    i53 = i53./mdet;
    i54 = i54./mdet;
    i55 = i55./mdet;

    
    Tdt  = -(i11.*R1+i12.*R2+i13.*R3+i14.*R4+i15.*R5);
    uxdt = -(i22.*R2+i23.*R3+i24.*R4+i25.*R5);
    uydt = -(i32.*R2+i33.*R3+i34.*R4+i35.*R5);
    P1dt = -(i42.*R2+i43.*R3+i44.*R4+i45.*R5);
    P2dt = -(i52.*R2+i53.*R3+i54.*R4+i55.*R5);


    % Tdt  = zeros(params.Nx, params.Ny);
    % uxdt = zeros(params.Nx, params.Ny);
    % uydt = zeros(params.Nx, params.Ny);
    % P1dt = zeros(params.Nx, params.Ny);
    % P2dt = zeros(params.Nx, params.Ny);

    % M = zeros(5, 5, params.Nx*params.Ny);
    % R = zeros(5, params.Nx*params.Ny);
    % M(1,1,:) = m11(:);
    % M(1,2,:) = m12(:);
    % M(1,3,:) = m13(:);
    % M(2,2,:) = m22(:);
    % M(2,3,:) = m23(:);
    % M(2,4,:) = m24(:);
    % M(2,5,:) = m25(:);
    % M(3,2,:) = m32(:);
    % M(3,3,:) = m33(:);
    % M(3,4,:) = m34(:);
    % M(3,5,:) = m35(:);
    % M(4,2,:) = m42(:);
    % M(4,3,:) = m43(:);
    % M(4,4,:) = m44(:);
    % M(5,2,:) = m52(:);
    % M(5,3,:) = m53(:);
    % M(5,4,:) = m55(:);
    % R(1,:) = R1(:);
    % R(2,:) = R2(:);
    % R(3,:) = R3(:);
    % R(4,:) = R4(:);
    % R(5,:) = R5(:);

    % sol = MultiSolver(M, R);



    timer4 = timer4 + toc;
    tic;




    udt = zeros(Nx.*Ny.*5, 1);

    udt(0*Nt+1 : 1*Nt) = Tdt(:);
    udt(1*Nt+1 : 2*Nt) = uxdt(:);
    udt(2*Nt+1 : 3*Nt) = uydt(:);
    udt(3*Nt+1 : 4*Nt) = P1dt(:);
    udt(4*Nt+1 : 5*Nt) = P2dt(:);

    timer5 = timer5 + toc();
    % udt(0*Nt+1 : 1*Nt) = gather(Tdt(:));
    % udt(1*Nt+1 : 2*Nt) = gather(uxdt(:));
    % udt(2*Nt+1 : 3*Nt) = gather(uydt(:));
    % udt(3*Nt+1 : 4*Nt) = gather(P1dt(:));
    % udt(4*Nt+1 : 5*Nt) = gather(P2dt(:));


end
